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A Parameter Choice Strategy for the Inversion of Multiple 

Observations 

Christian Gerhard^, Sergiy Pereverzyev JrJl, Pavlo TkachenkcH 
Abstract. 

In many geoscientific applications, multiple noisy observations of different origin need to be 
combined to improve the reconstruction of a common underlying quantity. This naturally 
leads to multi-parameter models for which adequate strategies are required to choose a set of 
’good’ parameters. In this study, we present a fairly general method for choosing such a set of 
parameters, provided that discrete direct, but maybe noisy, measurements of the underlying 
quantity are included in the observation data, and the inner product of the reconstruction 
space can be accurately estimated by the inner product of the discretization space. Then the 
proposed parameter choice method gives an accuracy that only by an absolute constant multi¬ 
plier differs from the noise level and the accuracy of the best approximant in the reconstruction 
and in the discretization spaces. 

Keywords. Parameter Choice, Multiple Observations, Spherical Approximation 


1 Introduction 

Satellite missions like CHAMP, GRACE, GOCE, or Swarm (e.g., [H El [TOl [13]) provide 
highly accurate data of the Earth’s gravity and magnetic field, e.g., by giving information 
on the first- or second-order radial derivative of the gravitational potential or measure¬ 
ments of the vectorial geomagnetic field, which, once certain iono- and magnetospheric 
contributions have been filtered out, can be expressed as the gradient of a harmonic 
potential. Drawing conclusions from such satellite measurements on the gravitational 
potential or the magnetic field at or near the Earth’s surface is a classical exponentially 
ill-posed problem (see, e.g., 1211191 [22]). Measurements at or near the Earth’s surface 
(which we simply denote as ground measurements), on the other hand, do not suffer from 
this ill-posedness but are typically only available in restricted regions (e.g., aeromagnetic 
surveying [23] ). Combining both data sets becomes necessary when aiming at local high 
resolution models that also take global trends into account. This is a classical setting 
for multiparameter modeling (e.g., |3l[I71[l8l[I2]) that involves the regularization of an 
ill-posed inverse problem (downward continuation of satellite data) and the weighting of 
the satellite data against the ground data. An exemplary situation that we also use for 
later numerical illustrations is the following: We have measurements fi of a harmonic 
potential u on a spherical satellite orbit D/j = {x € : |x| = i?} and measurements /2 
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of in a subregion C fir of the spherical Earth’s surface fir, r < R, i.e., 

Au = 0, in (1.1) 

u = / 2 , on Qr, (1.2) 

u = fi, on r^, (1.3) 

where = {x G : |x| > r}. The problem of approximating u in T^ is clearly 
overdetermined and spherical splines (e.g., la ED) or other localized basis functions 
(e.g., [I1I251I2S]) could be used to approximate u in T^ from knowledge of /i only 
(generally, we denote the restriction of u to T^ by v)). However, such methods are 
not always well-suited to capture global trends of u and they do not address situations 
where the noise level of /2 might be smaller than that of /i. Therefore, it is advisable to 
incorporate satellite data /2 as well. Eventually, based on different parameter settings or 
approximation methods, we assume to have a set of candidates {uk}k=i,2,...,N available 
for the approximation of u in 

In this paper, we aim at introducing a method that predicts a ’good’ candidate Uk* 
among the available {uk}k=i,2,...,N without requiring knowledge of the method by which 
each Uk has been obtained or which sort of noise is contained in the data. It is also 
not necessary to know the underlying models or the type of data that has lead to the 
construction of Uk- Apart from {uk}k=i, 2 ,...,N^ all that is required is a reference measure¬ 
ment / (in the example (in])-([L3l), this would be /i) of against which to compare the 
candidates Uk- In this sense, we are not dealing with a parameter choice strategy for an 
ill-posed problem (although the underlying models that determine u may be ill-posed) 
but rather with a general method of choosing a ’good’ approximant of among a set 
of available candidates (an extensive comparison of parameter choice methods for ill- 
posed problems can be found, e.g., in mm)- Opposed to aggregation methods (see, e.g., 
[3 ED]), where approximations from different data settings are superposed to obtain a 
hnal approximation, we assume in our method that this superposition has already taken 
place in one way or another during the construction of each Uk ■ An important constraint 
for our method, in order to obtain a suitable error estimate, is that the discrete refer¬ 
ence measurements / of need to be given in points that allow the definition of an 
inner product in the discretization space which coincides with the L^-inner product in a 
desired finite-dimensional function space (e.g., the spherical harmonics of degree smaller 
than some L). The numerical tests, however, show that our method also supplies good 
results if this condition is slightly violated. 

The structure of the paper is as follows: In Section [2l we introduce and investigate 
the parameter choice strategy mentioned above in more detail and put it into a mathe¬ 
matically rigorous context. In Section [S] we illustrate its performance for the problem 
(in])-((l3|). The approximations for this problem are obtained by a method described 
in m- Latter is also briefly recapitulated in Section [31 

2 The Parameter Choice Strategy 

Throughout this paper, we assume the following conditions to be satished: 
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(a) Let C be a subdomain of the sphere fir, where discrete direct measurements 

of the underlying quantity u are available. We assume to have M measurement 
values and a corresponding discretization operator D : L‘^[Tr) —>■ that maps a 

function € L‘^{Tr) to the correponding measurements Dtt^ G Furthermore, 
let be equipped with some inner product (•, •)]gM and the corresponding norm 

IMIkat • 

(b) The measurements of may be blurred by additive noise ^ = (^i, ■ ■ ■ ,Cm) £ 
and we assume, without loss of generality, that there is = / G L‘^{Tr) such that 


Du| = 


— Dut 


< e 


for some e > 0. 


(c) We assume that from somewhere, a set {ufc}fc=i 2 w approximations of on F^ 
is available sand that all these approximations belong to some hnite dimensional 
linear subspace V C L^(rr). 

(d) Finally, we assume that the discretization space is related to the reconstruction 
space V through the discretization operator D such that 


(9,9}L2(rr) = , for all g,g eV. (2.1) 


Example 2.1. Let V = Vl the space of spherical polynomials of the degree L. 
Under rather general assumptions on F^ one can find a system of knots ^ 

and positive weights such that 



M 

x)drr{x) = 

i=l 


for all g G V 2 L- 


Consider a discretization operator 


D 5 = {9{x^),9{x^), ■ ■ ■ ,g{x%)) G 


and the inner product 


M 

{y,y)RM ■■= 

i=l 


It is clear that for the just introduced reconstruction space V = Vl, discretization space 
M^, and discretization operator D the condition (|2.1D is satished. It is also clear that the 
measurements described by the operator D are just pointwise evaluations at the knots 
M' noise level of these measurements is controlled by the quantity 





M 

1=1 
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Now we are ready to describe our choice of a good approximation from a given family 
{uk}k=i,...,N- Let Uk,pt G {uk}k=i,...,N be such that 


- Uk, 


opt 


L^{Tr) 


mm 

k=l,2,...,N 


v) - Uk 


L^iTr) 


Of course, Uk^^t cannot be found without knowledge of . Therefore, in practice one 
cannot find the parameters corresponding to Uk^^^- 

We motivate our procedure with the observation that for any Uk, k = 1,2,N, it 
holds 


ll'^fc '^kopt\\L2(^p\ — sup {uk Uk„pt,a) 

aeL2(rr),||a|T2(p^)=l 

^ ~ L2(r^) ’ 

where the finite set is defined as follows 

An = \ a = ak,i = i ,———^ -, k,l = 1,2,..., n \ cV. 

[ W'^k J 

Then for any A; = 1, 2,..., and a G An the quantity 

{uk — '^kopti °)= {'^k,a)]^2(^Yr) ~ {'^kopti°) 

has only a part (uk^pt,'^) i, 2 (y ) that cannot be computed directly, because Uk^pt is un¬ 
known. On the other hand, this part can be approximated with the use of the available 
observations as follows 

= {Duk,p„Ba)^^ « (Du\,Da^^^ . 

Therefore, the values 

hk{a) = {uk, a) ^2^rr) “ 

and 


Hk = max \hk{a)\ 
aGAiv 

can be seen as surrogates for the values of (rtfc — Uk„pt,a)\\uk — Ukopt\\i2(^Y ) 
respectively. 

In the view of this it is natural to expect that the approximation Uk^ G {uk}k=i n 
defined by 

fe* : Hk, = min {Hk, k = 1,2,... , N} (2-2) 
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(2.3) 


is close to itfcopt. Indeed, 

\\uk, “ IIL2(rr) “ ~ '^kopti^k,,kopt) 

“ ((^feopt>®fc*,fcopt)i2(r^) - “fe*-fcopt)jgM) 

= /ifc* (Ofc*,fcopt) — hk^p^{(lkf,kopt) 

< + H}^ . < 2H]f .. 

- A/* ' f^opt - '^opt 


Furthermore, 


Hkopt = max 


= max 
cl^A^ 

<e + 


(^fcopt>a>r2m„^ - (Dn;,Do 


opt) “7^2 (r^) 

(o^fcopt - 

Du^ - Buk.pt 


Then the previous analysis together with the triangle inequality gives us the following 
statement. 


Theorem 2.2. Let us assume that conditions (a)~(d) hold true, i.e., we are given a 
finite family of approximations {uk}i^^i jv from a finite dimensional reconstruction 
space V C L‘^{Tr). Moreover, noisy direct discrete measurements Dtt^ G of the 
approximated quantity are available, and the reconstruction space V is related to the 
discretization space such that m is satisfied. Then for A:* chosen according to 

112.S^) we have 


- Uk, 


L^Tr) 


< 


- Uk 


opt 


L^Tr) 


+ 2 


Dtt^ — T)Uk, 


opt 


+ 26. 


Remark 2.3. Note that in the context of Example 12.II we can give also another bound 
for ||u^ — Ufc^||^ 2 (p )• Let G Vl be the spherical polynomial of the best C'(rp)- 
approximation, i.e.. 


Then 


^ O'hest 


= mm 
C(r,.) i;GVl 


u' — V 


CTr) 

— ( Dtt^, Da 


+ ( Dtt^ — Drtt, Da\ 

\ « /M*f 

= - 'wLea)^2(r ) “ 

+— Dut,Da\ + (uk 

\ V /rm \ ) / L2{Tr) 


<e + 


- Ukopt 


L^Tr) 


+ Cm,N 


- Ubest 


CiFr) 
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where 


M 


cm,n = Vvol(rr) + '^wf \a{xf)\ . 

i=l 

Furthermore, we observe that 

M / M \ / M \ 

i=l Vi=l / \i=l / 


= Vvol(f^. 


Thus, 


Hkopt < £ + 


v) - Uk, 


opt 


L^iFr] 


+ 2-\/vol(rr) 


- utst 


CiFr) 


and from (j2.3l) we get the following alternative bound: 


v} - Uk, 


< 3 


u' - Uk, 


opt 


L^iTr) 


+ 4y^vol(rr) 


^ ^best 


C(r,) 


+ 2s. 


Remark 2.4. For the estimates in Theorem 12.21 and Remark 12.31 the assumption (12.ip 
is crucial. Incorporating the worst-case error for non-exact quadrature rules in the 
reconstruction space V, estimates similar to those in Theorem 12.21 and Remark 12.31 conld 
be derived even if (12.11) is violated (compare, e.g., [15] for an overview on spherical 
quadratnre rules). However, such estimates would show a stronger and undesirable 
dependence on and Uk and are therefore omitted. In the numerical examples in the 
next section, we show that a ’slight’ violation of the condition 12.11 can still yield good 
results. 


3 Numerical Illustrations 

In this section, we illustrate the numerical performance of the previously described pa¬ 
rameter choice method on (inp-dLSI), i.e., we assume u to satisfy 

Au = 0, in 
u = / 2 , on Qr, 
u = fi, on F^. 

A set of approximations Uk, k = 1,... ,N, of on a spherical cap F^ = = {x G 

Hr : 1 — • (0,0,1)^ < p} of radius p G (0,2) around the North Pole (0,0, r)^ can be 

obtained by 

Uk{x)= ^k{x,y)f 2 {y)dfiRiy)+ ^k{x,y)fi{y)dTr{y), x G F^, (3.1) 

J JVr 
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where the kernels are given by 



with coefficients of the form ^^(n) = $^(n) — ‘h^(n) and N^. < M^. By 

{ynj]n=0,l,.. .;j=o,i,..., 2 n+i we mean a set of orthonormal spherical harmonics of degree 
n and order j. The coefficients <h^(n), ‘h^(n) are chosen by minimizing the following 
functional: 


Mk Nk 2 

|1 - 3.A(„)|2 + ^ |l _ 

n=0 n=0 

Nk 

+/3fc ^ I(^)I + ll^fc|L2(Q^\r^)- 

n=0 

The first two terms of the functional T measure the approximation property of the 
kernels 'I'fc (i.e., they measure how close they are to the Dirichlet kernel). The third 
term penalizes the error amplihcation due to the downward continuation of the satellite 
data on Qji while the fourth term penalizes the localization of outside the region 
Tr where ground data is available. The parameters ak, ak, Pk weigh these quantities 
against each other. For more details, on this approach of approximating u on T^, the 
reader is referred to m- Essentially, we are in the setting of Example 12.11 where the 
reconstruction space V = Vm is the space of all spherical polynomials up to degree 
M = max {Mfc, k = 1,..., A^}. 

The procedure for our numerical tests is as follows: 

(a) From the EGM2008 gravity potential model (cf. |2lj^ . we generated two sets of 
reference potentials u: 

(1) one up to spherical harmonic degree n = 30 (in order to allow many test 
runs in a short time) on a sphere Qr, R = 12,371km, and on a spherical cap 
Fr = Fr, r = 6,371km, with p = 1 (corresponding to a spherical radius of 
approximately 10,000km at the Earth’s surface), 

(2) another one up to spherical harmonic degree n = 130 (in order to have a 
more realistic scenario) on a sphere Qr, R = 7,071km, and on a spherical 
cap Fr = Fr, r = 6,371km, with p = 0.3 (corresponding to a spherical radius 
of approximately 5,000km at the Earth’s surface). 

^data accessed via http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html 
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(b) For both cases of part (a), we generate corresponding noisy measurements /i, 

/ 2 , where the noise levels ei = ||/i — ii||L2(rj.)/ll^llL2(rr) ground data and 

^2 = 11/2 “'u||L 2 (Q^)/||u||i 2 (Q^) of the satellite data are varied among 0.001, 0.1. 
The data on are in both cases computed on an equiangular grid according to 
[7] while the data on the spherical cap F^ are computed on a Gauss-Legendre grid 
according to [16] in order to guarantee polynomially exact quadrature rules up to 
spherical polynomial degree + n, where n = 30 in case (1) and n = 130 in case 
(2), which yields condition (|2.ip . 

(c) For the different input data from part (b), we compute approximations u^, k = 

1,..., A^, of on Fj. via the expression m- The index k of Uk indicates different 
choices of the parameters ak, dfc, /3k in the functional T from (l3|). Uk, ctk are 
varied in the interval [10^,10®] and /3k is varied in the interval [10“^,10^]. The 
truncation degrees of the series expansions of are fixed to Nk = Mk = 80 in 

case (1) while Nk = Mk = 150 in case (2). This way, we obtain N = 100 different 
approximations Uk for each of the two cases. 

(d) Among the approximations Uk, we use the procedure from Section [2] to choose 
a ’good’ approximation Uk*- Afterwards, we compare the relative approximation 
errors err^* = \\uk* — 'ii||L 2 (r^)/||u|| 2 , 2 (r^) of the parameter choice with the relative 
errors erropt = minfc=i,...,Ar \\uk - u||L 2 (r^)/||u||i 2 (r^) of the actually best Uk^^f 

The results of the tests are shown in Figures [1] and [2l Each figure shows the rela¬ 
tive errors err^* and erropt for every test run. Additionally, we plotted the maximum 
errors err^ax = maxk=i,...^N \\uk - u\\L2{rr)/\H\L2{rr-) ^nd the average errors errav = 
Jf Ylk=i N \Wk — '^‘lli, 2 (r^)/||tt|| 2 , 2 (r^) in order to illustrate the performance. It can be 
seen that the algorithm works particularly well for the setting £i = £2 and that the oracle 
error err^* is nearly identical with the minimum error erropt. The situation is different 
when £i S> £ 2 - The minimum error erropt is smaller than the noise level £ 1 . Thus, since 
our parameter choice strategy is based on comparing Uk to /i, we cannot expect that 
errfc* is as good as erropt. Yet, astonishingly enough, it seems that err^* is still slightly 
smaller than £1 for our test setting. 

In addition, we repeated the tests above with a reduced accuracy of the quadrature 
rule in order to illustrate the consequences if condition mi is not satisfied. More 
precisely, we did test runs for the following setting: 

(a’) We generated two sets of reference potentials u: 

(1) one up to spherical harmonic degree n = 30 on a sphere kin, R = 12, 371km, 
and on a spherical cap F^ = F^, r = 6,371km, with p = 1 (opposed to the 
previous tests, the potential is not based on the EGM2008 model but the 
Former coefficients are chosen randomly), 

(2) another one up to spherical harmonic degree n = 130 on a sphere fl/j, R = 
7,071km, and on a spherical cap F^ = F^, r = 6, 371km, with p = 0.3 (here, 
the potential is again based on the EGM2008 model). 
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Figure 1: Relative Errors for the low spherical harmonic degree tests (Situation (a)(1)) 
for £1 = £2 = 0.001 (left) and £i = 0.1, £2 = 0.001 (right; the dotted black line marks 
the noise level £1 = 0.1). 
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Figure 2: Relative Errors for the high spherical harmonic degree tests (Situation (a)(2)) 
for £1 = £2 = 0.001 (left) and £1 = 0.1, £2 = 0.001 (right; the dotted black line marks 
the noise level £1 = 0.1). 
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(b’) For both cases of part (a’), we generate corresponding noisy measurements /i, 
/2 with noise levels ei = £2 = 0.001. Again, the data on Qji are in both cases 
computed on an equiangular grid according to [7] while the data on the spherical 
cap Fr are computed on a Gauss-Legendre grid according to m- For case (1), we 
chose grids of two different sizes: one such that the polynomial exactness of the 
quadrature rule is of degree 100 and one such that polynomial exactness is of degree 
90 (remember that polynomial exactness up to degree + re = 80 + 30 = 110 is 
required in order to satisfy condition (|2.ip l. For case (2), we chose the size of the 
grids such that the polynomial exactness of the quadrature rule is of degree 130 
and of degree 80, respectively (remember that polynomial exactness up to degree 
Mk + re = 150 + 130 = 280 is required in order to satisfy condition (I2.ip l. 

(c’) For the different input data from part (b), we compute approximations Uk, k = 
1,... ,N, of on Fr via the expression (13.ip . The parameters otk are again 
varied in the interval [10^,10®] and f3k is varied in the interval [10“^,10^]. The 
truncation degrees of the series expansions of are fixed to = 80 

in case (1) while = 150 in case (2). 

The results are shown in Figures [3] and HI In the right plot of Figure [3] it becomes clear 
that a too large deviation of the required polynomial exactness can severely influence 
the parameter choice rule and render it essentially useless (the simple average of all ap¬ 
proximation errors is better than the error err^* of our algorithm). The left plot, on the 
other hand, shows that small deviations have hardly any influence. However, in order 
to illustrate this sensitive dependence on the polynomial exactness of the quadrature 
rule, we switched from the EGM2008 gravity potential model to potentials with Fourier 
coefficients that are generated randomly (i.e., in the mean, the Fourier coefficients are 
equally large at all spherical harmonic degrees). Figure H] shows that for a more re¬ 
alistic scenario like the EGM2008 model, the influence of the polynomial exactness of 
the quadrature rule is significantly smaller. In order to detect a severe failure of our 
algorithm, we had to decrease the polynomial exactness to degree 80 (opposed to degree 
280, which would guarantee the required condition (I2.ip l. This stability of the algorithm 
is due to the fact the the Fourier coefficients of the EGM2008 gravity potential show a 
fast decay for growing spherical harmonic degrees. The generally higher optimal errors 
eri’opt in Figures [31 [H compared to Figures [Il[2] have to be accounted to the influence of 
the decreased accuracy of the quadrature rule on the approximations Uk via (j3.ip but 
not to the parameter choice method presented in this paper. 

4 Conclusion 

We introduced a simple method to choose a ’good’ candidate Uk* among a set of approxi¬ 
mations {uk}k=i,...,N of and supplied some error estimates for u^* in relation to 
The numerical illustrations show its good performance and stability for applications, 
e.g., to the Earth’s gravity potential. 


10 


R=12742.4. 

0.04r 

0.035 [- 

0 . 03 |- 

0.025 [- 

0 . 02 [- 

0.015P 

0.011 

0.005p 

0*- 


r=6371.2, epsilon2=0.0010, epsilon^=0.0010, quadrature degree=100 



Q 

0 o 


R=12742.4, r=6371.2, epsilon2=0.0010, epsilon^=0.0010, quadrature degree=90 



Figure 3: Relative Errors for the low spherical harmonic degree tests (Situation (a’)(l)) 
for £i = £2 = 0.001 and a quadrature rule with polynomial exactness of degree 100 (left) 
and of degree 90 (right). 
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Figure 4: Relative Errors for the high spherical harmonic degree tests (Situation (a’)(2)) 
for £i = £2 = 0.001 and a quadrature rule with polynomial exactness of degree 130 (left) 
and of degree 80 (right). 
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